function [intensive,grossExtensive,selection] = calc_selection(Pdistans,Pdist_sh,lambda,D,model,interp_type,dh)

x0_ind = find(-Pdistans>=0,1);           %index of the first element, where gap is nonnegative
NNx = length(Pdistans);                       %number of gap-grids

posGapRange = x0_ind:NNx;                      
negGapRange = 1:x0_ind-1;

gapVec = -Pdistans;
posGap = -Pdistans(posGapRange)';
negGap = -Pdistans(negGapRange)';

%find F0
F_sh = cumsum(Pdist_sh);   %CDF of the after-shock distribution
F0 = interp1(gapVec,F_sh,0,interp_type);

%Approximate probabilities at gridpoints
posGap_mpoints = (posGap(1:end-1)+posGap(2:end))/2;
F_sh_posGap_mpoints = interp1(gapVec,F_sh,posGap_mpoints,interp_type);
HH_sh = ([F_sh_posGap_mpoints 1]-[F0 F_sh_posGap_mpoints])/(1-F0);       %probabilities at gridpoints

negGap_mpoints = (negGap(1:end-1)+negGap(2:end))/2;
F_sh_negGap_mpoints = interp1(gapVec,F_sh,negGap_mpoints,interp_type);
GG_sh = ([F_sh_negGap_mpoints F0]-[0 F_sh_negGap_mpoints])/F0;       %probabilities at gridpoints

LambdabarMinus = sum(lambda(posGapRange)'.*HH_sh);
LambdabarPlus = sum(lambda(negGapRange)'.*GG_sh);

intensiveMinus = LambdabarMinus;
intensivePlus  = LambdabarPlus;

inflMinus = sum(gapVec(posGapRange)'.*lambda(posGapRange)'.*HH_sh);
xbarMinus = sum(posGap.*HH_sh);
inflPlus = sum(gapVec(negGapRange)'.*lambda(negGapRange)'.*GG_sh);
xbarPlus = sum(negGap.*GG_sh);


switch model
 case 3 %Golosov-Lucas
  %Find lower and upper thresholds of price gaps

  %upper threshold xu
  yugrid = D(posGapRange)';
  xu = findinterp(posGap,yugrid,interp_type);   %upper threshold, where gain is 0

  %lower threshold xl
  ylgrid = D(negGapRange)';
  xl = findinterp(negGap,ylgrid,interp_type);   %lower threshold, where gain is 0

  f_sh = numderinterp1(gapVec,dh,gapVec,F_sh,interp_type);  %pdf of the after-shock distribution

  h_sh = f_sh/(1-F0);         %densities at gridpoints
  g_sh = f_sh/F0;

  h_sh_xu = interp1(gapVec,h_sh,xu,interp_type);  %density at the thresholds
  g_sh_xl = interp1(gapVec,g_sh,xl,interp_type);


  grossExtensiveMinus = xbarMinus*h_sh_xu;
  selectionMinus = (xu-xbarMinus)*h_sh_xu;

  grossExtensivePlus = -xbarPlus*g_sh_xl;
  selectionPlus = (xbarPlus-xl)*g_sh_xl;

 otherwise
  %derivative of the hazard
  %lambdaPrime = numderinterp1(gapVec,dh,gapVec,lambda,'pchip')';
  lambdaPrimeMinus = numderinterp1(gapVec(posGapRange),dh,gapVec(posGapRange),lambda(posGapRange),'pchip')';
  lambdaPrimebarMinus = sum(lambdaPrimeMinus.*HH_sh);
  grossExtensiveMinus = xbarMinus*lambdaPrimebarMinus;
  selectionMinus = sum(gapVec(posGapRange)'.*(lambdaPrimeMinus-lambdaPrimebarMinus).*HH_sh);
  
  lambdaPrimePlus = numderinterp1(gapVec(negGapRange),dh,gapVec(negGapRange),lambda(negGapRange),'pchip')';
  lambdaPrimebarPlus = sum(lambdaPrimePlus.*GG_sh);
  grossExtensivePlus = xbarPlus*lambdaPrimebarPlus;
  selectionPlus = sum(gapVec(negGapRange)'.*(lambdaPrimePlus-lambdaPrimebarPlus).*GG_sh);
end

intensive = F0*intensivePlus+(1-F0)*intensiveMinus;
grossExtensive = F0*grossExtensivePlus + (1-F0)*grossExtensiveMinus;
selection = F0*selectionPlus+(1-F0)*selectionMinus;

%test
%Fprime = @(x) numderinterp1(gapVec,dh,x,F_sh,interp_type);  %pdf of the after-shock distribution
%F = integral(Fprime,gapVec(1),gapVec(end));